-- Calculate an amplifier output pulse shape.
-- Requires XFun "fft".

decay = 1u     -- input exponential pulse decay time
bw1 = 200K     -- bandwidth of 1st amplifier stage
bw2 = 150K     -- bandwidth of 2nd amplifier stage

tmax = 8u      -- time scale
nsamp = 64     -- number of samples to use in transform

y(t) = exp(-t/decay) -- input pulse

Xmin=0; Xmax=tmax/1u; t = X*1u; Xlabel="sec"
plotline y(t)

newaxis
dt=tmax/(nsamp-1)
ys[i]={y((i-1)*dt),0} -- sample the input to get a complex time series
fft(ys[:nsamp]):;     -- transform the input signal to freq domain
fmax=.5/dt;  fmax:3.94e6; -- highest freq used in model

H(f,bw) = Cdiv({1,0},{1,f/bw}) -- freq response of a single RC lowpass
freq(i) = (i-1)*fmax/(nsamp/2+1)  
amp(bw)[i] = H(freq(i),bw)  -- resp at discrete freqs of the transform

stage1:=convolve(transform,  -- convolve transform of input pulse
                 amp(bw1)):; -- with 1st amplifier response
invfft(stage1):;             -- invert to get time series
plotline realpart(transform) -- show 1st stage output pulse shape

invfft(convolve(stage1,      -- convolve 1st stage output freqs
                amp(bw2))):; -- with 2nd amplifier response
plotline realpart(transform) -- show 2nd stage output pulse shape


------------ utility routines for complex numbers ----------------
Cmult(A,B) = {A[1]*B[1]-A[2]*B[2],A[2]*B[1]+A[1]*B[2]}
Cdiv(A,B)={(B[1]*A[1]+B[2]*A[2])/(B[1]^2+B[2]^2),
         (B[1]*A[2]-B[2]*A[1])/(B[1]^2+B[2]^2)}
convolve(A,B)[i,j] = Cmult(A[i],B[i])[j] dim[count(A),2]
realpart(A)[i] = A[i,1] dim[count(A)]
